SAIC-01 /8020 : APPAT-287 
November 26, 2001 


"Global magnetohydrodynamic 
Modeling of the Solar Corona" 

NASW-98030 
Final Report, 

Sept. 30, 1998 - Sept 29, 2001 


NASA CONTRACT: NASW-98030 


PRINCIPAL INVESTIGATOR: 

JON A. LINKER 

SCIENCE APPLICATIONS INTERNATIONAL CORPORATION 

10260 Campus point Drive 
San Diego, CA 92121-1578 


SAIC-01/8020:APPAT-287 
November 26, 2001 

GLOBAL MAGNETOHYDRODYNAMIC MODELING OF THE SOLAR CORONA 

FINAL REPORT: 

9/30/98 - 9/29/01 


1. Introduction 

The solar corona, the hot, tenuous outer atmosphere of the Sun, exhibits many 
fascinating phenomena on a wide range of scales. One of the ways that the Sun can 
affect us here at Earth is through the large-scale structure of the corona and the 
dynamical phenomena associated it, as it is the corona that extends outward as the 
solar wind and encounters the Earth's magnetosphere. The goal of our research 
sponsored by NASA's Supporting Research and Technology Program in Solar 
Physics is to develop increasingly realistic models of the large-scale solar corona, so 
that we can understand the underlying properties of the coronal magnetic field that 
lead to the observed structure and evolution of the corona. In the following 
sections we describe the work performed under this contract. 


2. Computational Modeling of the Corona and Inner Heliosphere 


Support from this contract has allowed us to develop sophisticated 
magnetohydrodynamic (MHD) computations of the solar corona and inner 
heliosphere. A key aspect of the computations are that observed photospheric 
magnetic fields are incorporated into the boundary conditions, allowing the 
performance of computations that can be compared with specific observations. To 
compute self-consistent three-dimensional MHD solutions for the large-scale 
corona, we solve the following equations in spherical coordinates: 
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+ V-(pv) = (y - 1)( - pV-v + 5) , (6) 

S = - Vq - n e n p Q(T) + H c h + Hd + D , (7) 

where B is the magnetic field, J is the current density, E is the electric field, p, v, p, 
and T are the plasma mass density, velocity, pressure, and temperature, and the 
wave pressure p w represents the acceleration due to Alfven waves. The 
gravitational acceleration is g, y = 5/3 is the ratio of specific heats, t) is the resistivity, 
v is the viscosity, H c h is the coronal heating source, D is the Alfven wave dissipation 
term, n e and rip are the electron and proton density, and Q(T) is the radiation loss 
function (Rosner et al. 1978). The term Hd = r)J 2 + vVv:Vv represents heating due to 
viscous and resistive dissipation. In the collisional regime (below ~ 10R S ), the heat 
flux is q = -X| |bb-VT, where b is the unit vector along B, and K| j =9 x 10‘ 7 T 5 / 2 is 
the Spitzer value of the parallel thermal conductivity. In the collisionless regime 
(beyond ~ 10R S ), the heat flux is given by q = an e kT\, where a is a parameter 
(Hollweg 1978). Since it is presently not known in detail what heats the solar 
corona, the coronal heating source H c h is a parameterized function (see Lionello, 
Linker, and Mikic, 2001). Our model can also incorporate the acceleration of the 
solar wind by Alfven waves using a WKB approximation (Hollweg, 1978) as 
described by Mikic et al. (1999). Note that the simplified polytropic model is obtained 
by setting S = 0 in Eq. (6), p w = 0 in Eq. (5), and y = 1.05. 

The methods we use to solve equations (1-7), including the boundary conditions, 
have been described previously (Mikic and Linker 1994, 1996; Linker et al., 1996; 
Linker and Mikic, 1997; Mikic et al., 1999; Lionello, Mikic, and Linker, 1999). The 
principal observational input to our MHD model is the measured line-of-sight 
photospheric magnetic field. Typically, we have used synoptic magnetic field maps, 
which are generated from daily measurements of the magnetic field on the visible 
solar disk in the photosphere and are available from a number of observatories 
(including Wilcox Solar Observatory, Mt. Wilson Solar Observatory, the National 
Solar Observatory (NSO) at Kitt Peak and the MDI instrument aboard SOHO). 

The radial component of the field ( B r0 ) is deduced from the measured line-of- 
sight component (Wang and Sheeley 1992). In addition to B r0 we must specify the 
density ( p r0 ) and temperature (T r0 ) as boundary conditions at the lower boundary. A 
potential magnetic field consistent with the specified B^and a transonic Parker solar 
wind solution consistent with the specified p w and T r0 are used as initial conditions. 
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and Equations (1-7) are integrated in time until a steady-state is reached. The 
plasma flow velocity parallel to the magnetic field at the lower boundary is not 
specified but is obtained using the MHD characteristic equations. Characteristic 
equations are also used to calculate quantities at the outer boundary (typically placed 
at 30R S ), where the flow is supersonic and super-Alfvenic. The solution obtained 
provides a 3D description of the solar corona, under the assumption that the large- 
scale coronal magnetic field is not changing significantly for the time period of 
interest. 

3. Comparisons of the Polytropic Model with Coronal and Heliospheric Data 

When S is set to zero, equations (1-7) yield polytropic solutions. These 
solutions have the advantage that relatively simple models can match many of the 
properties of the corona; however, values of y close to 1 (y = 1.05 for the results 
shown here) are necessary to produce plasma profiles that are similar to coronal 
observations (Parker, 1963). Using the plasma density from these computations, we 
have computed the polarization brightness (pB) as would be observed from Earth or 
from spacecraft and developed simulated eclipse and coronagraph images that can 
be compared with corresponding measurements. We have performed extensive 
comparison of the polytropic model with eclipse observations (Mikic and Linker 
1996; Linker et al., 1996; Linker and Mikic, 1997; Mikic et al., 1999), and posted 
predictions and comparisons of coronal structure on our web site 
(http: / /haven.saic.com / corona / modeling.html) for the last 5 total solar eclipses. 

We have also compared our results with Mauna Loa Coronameter images, LASCO 
images, coronal holes as deduced from EIT images, and NSO Kitt Peak coronal hole 
maps (Linker et al., 1999ab). Linker et al. (1999b), Gibson et al. (1999), and Breen et al. 

(1999) also describe comparisons of our MHD model with observations during 
Whole Sun Month (WSM; August 10 - September 8, 1996) when a wide range of 
coronal and heliospheric data was available for coordinated study. 

The comparisons discussed above were primarily for times near solar 
minimum. To extend our technique to compute coronal solutions during solar 
maximum, two difficulties present themselves. The first problem is that the 
photospheric magnetic field is changing more rapidly, so synoptic magnetic maps 
may poorly represent the photospheric magnetic field. The second difficulty is that 
much higher resolution is required in an MHD computation to capture the highly 
structured, complex features of the solar maximum corona. Figure 1(a) shows 
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magnetic field lines for an MHD computation for Carrington rotation 1951 (June 24-- 
July 21, 1999). This calculation was performed in preparation for our prediction of 
coronal structure during the August 11, 1999 solar eclipse, on a 111 x 141 x 128 
nonuniform ( r,Q,(p ) mesh. Note the presence of helmet streamers (regions of closed 
magnetic field) at nearly all latitudes and the complex nature of the photospheric 
magnetic field (shown on the solar surface). The predicted radial plasma flow 
velocity is contoured in the plane of the sky for this view in Figure 1(b). We see that 
the slow flow regions occur at much higher latitudes in comparison to solar 
minimum (for example, see Plate 2 of Linker et al., 1999b). Figure 2 shows a 
comparison of simulated pB images with data from the Mauna Loa MKIII 
coronameter. Despite the difficulties of computing solar maximum solutions, the 
MHD model appears to have captured many features of the overall structure of the 
corona during this time period. However, the Mauna Loa images clearly show 
more fine structure than is present in the MHD computation, and we also found 
that the comparisons became worse as the difference between the time of the 
photospheric measurements and the time of the white-light observations increased. 
Mikic et al. (2000) describes this calculation in more detail, including our 
comparison with eclipse images. 

We have also extended our MHD calculations of the solar corona out into the 
inner heliosphere. We have found that the most efficient way to accomplish these 
computations is to split the calculation into a coronal solution (typically 1-30 R^ and 
a heliospheric solution (30 R s to 5 A.U.) Even with the use of implicit methods, we 
find that the inner solution requires a time step that is prohibitively small for the 
heliospheric part of the calculation. A more economical approach is to compute a 
separate heliospheric solution with a much larger time step. As long as the location 
of the interface between the calculations is beyond the Alfven and sonic points, it is 
relatively straightforward to use values from the inner coronal solution as the 
boundary values for the heliospheric calculation (all the MHD characteristics point 
outward into the heliosphere, so no information propagates back upstream). W e 
use the magnetic field computed from the coronal solution to specify the boundary 
condition for the outer solution. At the present time, because the velocity from our 
polytropic solution does not yield the observed contrast between fast and slow solar 
wind, we prescribe the velocity at 30 R s using the topology of the magnetic field (fast 
flow is assumed to originate from deep within coronal holes, while slow wind is 
assumed to come from the boundaries of coronal holes). This relatively simple 
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Figure 1. (a) Magnetic field lines from a computation of coronal structure near solar maximum 
(Carrington rotation 1951; June 24 July 21, 1999). Contoured on the solar surface is the 
photospheric magnetic field deduced from synoptic magnetic field measurements. In contrast to solar 
minimum, helmet streamers can be seen at nearly all latitudes, (b) Color contours of the plasma 
radial velocity from the model, in the plane of the sky of Figure la. Note that the regions of stagnant 
and slow flow extend to very high latitudes. 


June 25, 1999 June 28, 1999 July 2, 1999 July 6, 1999 



Figure 2. Comparison of results from the MHD model shown in Figure 1 with white light 
observations. The top panels show the polarization brightness (pB) calculated using the plasma 
density predicted by the MHD model, and the bottom panels show the observed pB from the Mauna 
Loa MKIII coronameter. The observed images in general show more structure than is present in the 
computations, and there are individual discrepancies, but the MHD model appears to have captured 
the overall structure of the corona, during this interval. 
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Figure 3. 3D MHD model of the heliosphere during Whole Sun Month (August 10 September 8,1996). 
(a) Meridional (plane of constant longitude) view of the radial velocity, scaled plasma number density, 
and scaled plasma pressure, (b) Comparison of Ulysses observations of solar wind speed, proton 
number density, and proton temperature with the MHD model. 
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prescription leads to the specification of a complex speed profile at the inner 
boundary of the heliospheric MHD calculation. (Eventually, the thermodynamic 
solutions described in section 2.3 will allow us to specify the plasma parameters 
directly from the MHD solution as well.) 

Figure 3(a) shows meridional cuts (planes of longitude) for a calculation of 
the inner heliosphere performed for the Whole Sun Month time period. Color 
contours of the radial velocity, scaled plasma number density (to account for the 
rapid decrease of density with radial distance), and scaled plasma pressure are shown 
along with the position of the heliospheric current sheet (white line). Note the 
steepening of the density and pressure that occurs when fast wind catches up with 
slow wind. During WSM, Ulysses was located at 4.25 AU from the Sun, at a 
heliographic latitude of 28°N. Figure 3(b) compares Ulysses observations (blue) of 
bulk flow speed, proton number density, and proton temperature with the 
simulation results (red), which were obtained by flying the Ulysses trajectory 
through the model. Overall, the simulation reproduces the large-scale features 
observed during this time period. Riley et al. (2001ab) describe these and other data 
comparisons in more detail. Posner et al. (1999,2001) describes additional 
comparisons of our model interplanetary in situ data, and Breen et al. (1999) 
discusses comparisons of the model with interplanetary scintillation (IPS) 
measurements. 

We have used our 3D MHD model of the corona and inner heliosphere to 
explore evolution of the heliosphere over the course of the solar cycle. To illustrate 
the evolution of the large-scale structure of the HCS during the course of the solar 
cycle, we present a summary of several solar parameters, measured over a period of 
-2-1/2 cycles, together with a selection of simulation results in Figure 4. The central 
panel includes data from cycles 21, 22, and the ascending phase of 23 and contains: 

(1) The average computed tilt angle (black) of the HCS as derived from source- 
surface calculations using photospheric measurements from the Wilcox Solar 
Observatory by T. Hoeksema (http: / / quake.stanford.edu / ~wso / Tilts.html ) and 
smoothed over 3 Carrington rotations; and (2) Monthly (yearly)- averaged values of 
Sunspot number in red (blue). Note the significant correlation between the WSO- 
derived tilt angle and the smoothed sunspot number: Both show a faster rise during 
the ascending phase of the solar cycle and a slower decay during the descending 
phase. 
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Figure 4. Evolution of several solar parameters during cycles 21, 22 , and the 
ascending phase of 23 , with emphasis on the evolution of the HCS during solar 
cycle 22. The lower-central panel shows the monthly (yearly )-averaged values of 
Sunspot number in red (blue). The upper-central panel shows the m=0 
azimuthally-symmetric part of the radial component of the magnetic field, as 
inferred from Kitt Peak synoptic maps, for each Carrington rotation. Blue 
indicates inward polarity and red indicates outward polarity. The HCSs as 
computed by the MHD model (out to 5 AU) for eleven Carrington rotations, 
covering mid-1986 to mid-1996, are shown above and below the central panels. 
Time runs from the top left to right, and then bottom left to right. 
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The HCS for eleven Carrington rotations, covering mid-1986 to mid-1996, are 
shown above and below the central panels. Time runs from the top left to right, and 
then bottom left to right and each simulation is separated from the next by ~13 
Carrington rotations. These isosurfaces are in qualitative agreement with the WSO- 
derived tilt angle profile. In particular, the rapid growth of the HCS to high 
heliographic latitudes during the ascending phase and the slower decline during the 
descending phase. A more detailed inspection of the isosurfaces reveals several 
noteworthy features. First, surrounding solar minimum, the HCS is better described 
as a flat surface with one or more folds in it, in contrast to the sinusoidal picture that 
is generated by considering the interplanetary extension of a tilted dipole. Second, 
folds of the HCS are typically asymmetric with respect to heliocentric distance: A 
fold rises more sharply on the inner radial side and falls more slowly on the outer 
radial slide. This is a natural consequence of the dynamic interaction of the 
surrounding streams and is particularly effective near solar minimum. Adopting a 
simplified picture of slow solar wind flow being organized about the HCS, and faster 
flow elsewhere, this asymmetry can be understood by the "stretching" of the HCS on 
the outer portion of the fold, where a rarefaction region is developing and slower 
flow is being accelerated into it, and a "squashing" of the HCS on the inner portion 
of the fold due to the formation of a compression region. 

4. Improvements to the Energy Equation 

While the results from the poly tropic model are encouraging, we know from 
detailed comparisons that the model is not sophisticated enough to yield the strong 
temperature variations that are observed in corona. The plasma density from the 
polytropic model can appear similar to the observations qualitatively, but the 
density contrast between coronal holes and the streamer belt is not accurate 
quantitatively, nor is the solar wind velocity far from the Sun. These difficulties 
motivated us to formulate a more realistic energy equation (including the terms in 
S shown in eq. (7)), similar to what had been performed in previous one- 
dimensional studies (Withbroe, 1988). The inclusion of these energy transport 
processes into our MHD model allow us to place the lower boundary in the upper 
chromosphere (T = 20,000 K), and still compute a solution of the large-scale solar 
corona. We refer to this as our chromospheric MHD model. Figure 5 shows 
magnetic field lines and the plasma temperature from a two-dimensional 
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Figure 5. The plasma temperature (T) in the solar corona from an MHD simulation that includes the upper 
chromosphere and transition region. Also shown are magnetic field lines. Blue shows the lowest temperatures 
and red the highest. T varies from less than 20,000 K in the upper chromosphere to more than 2,000,000 K in 
the corona. 
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Figure 6. The plasma number density as a function of radial distance for the polytropic and chromospheric 
MHD models. Note that the chromospheric model incorporates the steep temperature gradient in the lower 
transition region, and that a much larger contrast between the polar and equatorial density occurs in this 
model.The inset images show the polarization brighness computed for each model. Note the larger contrast 
in brightness between the streamer and the polar coronal holes for the chromospheric MHD model. 



Figure 7. Comparison of EIT images of an equatorial coronal hole on August 27, 1996 with MHD models, 
(a) Open field regions (black) from the polytropic MHD Model.(b) EIT 195A image showing the coronal 
hole, (c) EIT 284A image (d) Simulated EIT 284A image using the thermodynamic MHD model. 
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(azimuthally symmetric) calculation as an illustration. The temperature can be seen 
to vary on a wide range of scales, but is still captured by the model. To capture these 
sharp density and temperature gradients, the smallest grid cells in our calculation 
were 70km, or 0.0001 R s (a nonuniform 200 x 300 r,Q mesh was used). These results 
are described in more detail by Lionello et al. (2001). 

To illustrate the quantitative differences in solutions with the improved 
model versus the polytropic model. Figure 6 shows radial profiles of the plasma 
density at the equator and poles for the two different models. In the chromospheric 
model the equatorial (streamer density) is nearly an order of magnitude greater than 
the polar (coronal hole density), similar to observations. The polytropic model 
shows a much weaker density contrast between the equator and pole. The insets 
show simulated polarization brightness images for these two models. Note that the 
polar coronal holes are quite dark in the chromospheric MHD model (as is 
observed). The improved thermodynamic description in our MHD model makes 
modeling disk emission possible, just as we have previously done for polarization 
brightness (pB). The more realistic temperature obtained from the solution can be 
used to predict the abundance of the coronal iron species (using CHIANTI, Dere et 
al., 1997) and produce "simulated" EIT images. To illustrate the idea, we show a 
comparison we performed for the Whole Sun Month period (Linker et al., 1999b). 

Open field regions predicted by the polytropic MHD model for August 27, 1996 (Fig. 

7a) are compared with EIT images showing an equatorial coronal hole on that day 
(Fig. 7b and c). For our first 3D computation with our improved thermodynamic 
model we recalculated the WSM case, and developed a "simulated" EIT image 
(Figure 7d). In this preliminary case we did not include the upper chromosphere, 
and there is insufficient resolution to reproduce the fine structure seen in EIT, but 
the model does show low emission in the vicinity of the coronal hole. TRACE 
emission lines are similar to EIT, and we can use a similar procedure to develop 
simulated Yohkoh images. In the future, we plan to use simulated emission images 
to test our models more extensively. 

5. Modeling Prominence Formation within a Helmet Streamer 

One of our motivations for including the chromosphere is to allow the 
possibility of realistically modeling prominence formation. As a first step in this 
type of calculation, we have used a 2.5D axisymmetric MHD model to self- 
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Figure 8. The evolution of the plasma temperature and magnetic field in a 
thermodynamic MHD model. The plasma temperature is depicted in color, and 
projections of the magnetic field lines are overlaid on the temperature. A closeup of the 
lower part of the simulation domain near the equator is shown. (a) The helmet 
streamer configuration at the end of the shearing phase, (b) The streamer after the 
magnetic flux near the neutral line has been reduced 3.75% of its initial value. A low- 
lying filament structure has formed and is just discernible near the lower boundary, 
(c) The filament is now at a height of 140,000 km and is moving upward slowly. The 
enhanced density can be seen to lie near the bottom of the detached flux surfaces. Flux 
reduction (at 11.25%) is beyond the critical threshold for eruption (see text), (d) the 
eruptive phase has started, and dense material is carried into the outer corona, shown 
in (e) and (f). 
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consistently describe the formation of a stable prominence that supports cool dense 
material in the lower corona. Reducing the magnetic flux along the neutral line of a 
sheared coronal arcade creates a magnetic field configuration with a flux-rope 
topology. Formation of the flux lifts dense chromospheric material into the corona. 
This prominence-like structure sits at the base of a helmet streamer structure. The 
dense material is supported against gravity in the dips of the magnetic field lines in 
the flux rope. Further reduction in magnetic flux leads to an eruption of the 
prominence, ejecting material into the solar wind. Figure 8 shows the plasma 
temperature during the formation and ejection of the prominence. A paper 
describing these results is in press in the Journal of Geophysical Research (Linker et 
al., 2001). 

6. Summary 

Support from the NASA SR&T contract "Global Magnetohydrodynamic 
Modeling of the Solar Corona" has allowed us to greatly advance the state of 
computational modeling of the solar corona. We have extensively compared our 
computational models of the corona with eclipse and ground-based coronameter 
observations (Mikic et al., 1999, 2000; Linker et al., 1999a), observations from the 
Large Angle Spectroscopic Coronagraph (LASCO) and the Extreme ultraviolet 
Imaging Telescope (EIT) aboard the Solar and Heliospheric Observatory (Linker et 
al., 1999b; Gibson et al., 1999) and interplanetary IPS and in situ measurements 
(Breen et al., 1999; Posner et al., 1999, 2001; Riley et al., 2001ab). These comparisons 
have shown that MHD models can capture many important features of coronal and 
heliospheric structure, including the location and shape and of the coronal streamer 
belt, the distribution of coronal holes, the location of fast and slow solar wind 
streams, and the variation of coronal and heliospheric structure as a function of 
solar cycle. Inaccuracies of the model in computing coronal temperatures and solar 
wind velocities motivated us to improve the energy equation in the model and to 
simulate processes deeper in the solar atmosphere. We have extended the physics 
in our MHD model to include parallel thermal conduction, coronal heating, and 
radiation losses in the energy equation (Mikic et al., 1999). These computations can 
now simulate the top of the chromosphere (at 20,000 K) through the transition 
region as well as the corona and inner heliosphere (Lionello et al., 2001). We have 
applied this capability to modeling of prominence formation within a helmet 
streamer (Linker et al., 2001). A list of our publications follows in the next section. 
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